logo_datactivist           logo_pleiade


Objectif : Exploration des données et modélisations

Date début de l’analyse : 17 mai 2022

Date de la dernière modification : 24 mai 2022

# Packages
packages = c("tidyverse", "glue", "kableExtra", "knitr", "plotly", "hrbrthemes", "gridExtra", "stats", "arrow", "lme4", "tidymodels", "broom.mixed", "multilevelmod", "dotwhisker", "rAmCharts", "DescTools")

## Installation des packages si besoin et chargement des librairies
package.check <- lapply(packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)}})

# Import BSO complet (2013-2020)
data <- read_csv("../data/4.process/managing-variables_2nd-part/whole_bso.csv", col_types = cols(author_position = "c"), na = c("", "NA"))

# Création de nouvelles variables : couleur de journal finale et APC fiables
data <- data %>% mutate(oa_color_finale = case_when(oa_color_article_BSO == "diamond" ~ "diamond",
                                                    oa_color_openalex == "gold" ~ "gold",
                                                    oa_color_openalex == "hybrid" ~ "hybridOA",
                                                    TRUE ~ NA_character_), #qs Maurits : catégorie "other" ou tout en NA ? Sur quelle variable se baser alors ?
                        amount_apc_EUR_trusted = case_when(apc_source == "openAPC_estimation_publisher" | apc_source == "openAPC_estimation_publisher_year" ~ NA_real_,
                                                           TRUE ~ amount_apc_EUR))

I - Analyse exploratoire


Y moyen par année et par variable catégorielle


Y1 = amount_apc_EUR

Discipline

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>%  group_by(year, bso_classification) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup()

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = bso_classification, colour = bso_classification)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par discipline et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Tier

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>% group_by(year, tier) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup() %>% mutate(tier = as.character(tier))

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = tier, colour = tier)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par tier et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


French CA

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>% 
    group_by(year, is_french_CA_bso_wos_openalex_single_lang) %>% 
    summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup() %>% 
    mutate(is_french_CA = as.character(is_french_CA_bso_wos_openalex_single_lang))

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_french_CA, colour = is_french_CA)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si le CA est français") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Couleur OA

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>% group_by(year, oa_color_finale) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup() %>% mutate(oa_color_finale = as.character(oa_color_finale))

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = oa_color_finale, colour = oa_color_finale)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par couleur OA et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Covered by Couperin

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR)) %>% group_by(year, is_covered_by_couperin) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR, na.rm = T)) %>% ungroup() %>% mutate(is_covered_by_couperin = as.character(is_covered_by_couperin))

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_covered_by_couperin, colour = is_covered_by_couperin)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si l'article est couvert par les accords Couperin") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Y2 = amount_apc_EUR_trusted

Discipline

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>%  group_by(year, bso_classification) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = bso_classification, colour = bso_classification)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par discipline et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :

  • doaj ;
  • OpenAPC ;
  • openAPC_estimation_issn_year ;
  • openAPC_estimation_issn.

Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.

# Sélection de variables pour les modèles
bso_mvp <- data %>% 
    select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>% 
    mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.factor(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA_bso_wos_openalex_single_lang = as.factor(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.factor(is_covered_by_couperin))

# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>% 
    filter(!is.na(bso_classification),
           !is.na(oa_color_finale),
           !is.na(is_french_CA_bso_wos_openalex_single_lang),
           !is.na(amount_apc_EUR_trusted),
           !is.na(is_covered_by_couperin)) %>% 
    filter(apc_has_been_paid == 1, amount_apc_EUR_trusted > 0)

À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :

  • sélection : sélection des variables qui entreront dans les modélisations ;
  • filtres : premier filtre pour ne garder que les articles pour lesquels toutes les variables à intégrer aux modèles ont une valeur connue (not NA), puis second filtre sur les articles pour lesquels des APC ont été payés (apc_has_been_paid == 1, amount_apc_EUR_trusted > 0).

Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.


Tier

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, tier) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(tier = as.character(tier))

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = tier, colour = tier)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par tier et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :

  • doaj ;
  • OpenAPC ;
  • openAPC_estimation_issn_year ;
  • openAPC_estimation_issn.

Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.

# Sélection de variables pour les modèles
bso_mvp <- data %>% 
    select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>% 
    mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.factor(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA_bso_wos_openalex_single_lang = as.factor(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.factor(is_covered_by_couperin))

# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>% 
    filter(!is.na(bso_classification),
           !is.na(oa_color_finale),
           !is.na(is_french_CA_bso_wos_openalex_single_lang),
           !is.na(amount_apc_EUR_trusted),
           !is.na(is_covered_by_couperin)) %>% 
    filter(apc_has_been_paid == 1, amount_apc_EUR_trusted > 0)

À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :

  • sélection : sélection des variables qui entreront dans les modélisations ;
  • filtres : premier filtre pour ne garder que les articles pour lesquels toutes les variables à intégrer aux modèles ont une valeur connue (not NA), puis second filtre sur les articles pour lesquels des APC ont été payés (apc_has_been_paid == 1, amount_apc_EUR_trusted > 0).

Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.


French CA

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% 
    group_by(year, is_french_CA_bso_wos_openalex_single_lang) %>% 
    summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% 
    mutate(is_french_CA = as.character(is_french_CA_bso_wos_openalex_single_lang))

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_french_CA, colour = is_french_CA)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si le CA est français") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :

  • doaj ;
  • OpenAPC ;
  • openAPC_estimation_issn_year ;
  • openAPC_estimation_issn.

Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.

# Sélection de variables pour les modèles
bso_mvp <- data %>% 
    select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>% 
    mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.factor(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA_bso_wos_openalex_single_lang = as.factor(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.factor(is_covered_by_couperin))

# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>% 
    filter(!is.na(bso_classification),
           !is.na(oa_color_finale),
           !is.na(is_french_CA_bso_wos_openalex_single_lang),
           !is.na(amount_apc_EUR_trusted),
           !is.na(is_covered_by_couperin)) %>% 
    filter(apc_has_been_paid == 1, amount_apc_EUR_trusted > 0)

À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :

  • sélection : sélection des variables qui entreront dans les modélisations ;
  • filtres : premier filtre pour ne garder que les articles pour lesquels toutes les variables à intégrer aux modèles ont une valeur connue (not NA), puis second filtre sur les articles pour lesquels des APC ont été payés (apc_has_been_paid == 1, amount_apc_EUR_trusted > 0).

Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.


Couleur OA

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, oa_color_finale) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(oa_color_finale = as.character(oa_color_finale))

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = oa_color_finale, colour = oa_color_finale)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par couleur OA et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :

  • doaj ;
  • OpenAPC ;
  • openAPC_estimation_issn_year ;
  • openAPC_estimation_issn.

Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.

# Sélection de variables pour les modèles
bso_mvp <- data %>% 
    select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>% 
    mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.factor(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA_bso_wos_openalex_single_lang = as.factor(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.factor(is_covered_by_couperin))

# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>% 
    filter(!is.na(bso_classification),
           !is.na(oa_color_finale),
           !is.na(is_french_CA_bso_wos_openalex_single_lang),
           !is.na(amount_apc_EUR_trusted),
           !is.na(is_covered_by_couperin)) %>% 
    filter(apc_has_been_paid == 1, amount_apc_EUR_trusted > 0)

À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :

  • sélection : sélection des variables qui entreront dans les modélisations ;
  • filtres : premier filtre pour ne garder que les articles pour lesquels toutes les variables à intégrer aux modèles ont une valeur connue (not NA), puis second filtre sur les articles pour lesquels des APC ont été payés (apc_has_been_paid == 1, amount_apc_EUR_trusted > 0).

Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.


Covered by Couperin

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, is_covered_by_couperin) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(is_covered_by_couperin = as.character(is_covered_by_couperin))

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_covered_by_couperin, colour = is_covered_by_couperin)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si l'article est couvert par les accords Couperin") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Au vu des graphiques présentés ci-dessus, nous décidons de ne modéliser que les articles pour lesquels le montant des APC est fiable, c’est-à-dire provenant de l’une des sources suivantes :

  • doaj ;
  • OpenAPC ;
  • openAPC_estimation_issn_year ;
  • openAPC_estimation_issn.

Parmi les 6 sources d’APC recensées par le BSO, nous écartons ainsi “openAPC_estimation_publisher” et “openAPC_estimation_publisher_year”.

# Sélection de variables pour les modèles
bso_mvp <- data %>% 
    select(year, bso_classification, apc_has_been_paid, tier, oa_color_finale, is_french_CA_bso_wos_openalex_single_lang, amount_apc_EUR_trusted, is_covered_by_couperin) %>% 
    mutate(bso_classification = as.factor(bso_classification), apc_has_been_paid = as.factor(apc_has_been_paid), tier = as.factor(tier), oa_color_finale = as.factor(oa_color_finale), is_french_CA_bso_wos_openalex_single_lang = as.factor(is_french_CA_bso_wos_openalex_single_lang), is_covered_by_couperin = as.factor(is_covered_by_couperin))

# Préparation du jeu de données prêt à modéliser
bso_pop <- bso_mvp %>% 
    filter(!is.na(bso_classification),
           !is.na(oa_color_finale),
           !is.na(is_french_CA_bso_wos_openalex_single_lang),
           !is.na(amount_apc_EUR_trusted),
           !is.na(is_covered_by_couperin)) %>% 
    filter(apc_has_been_paid == 1, amount_apc_EUR_trusted > 0)

À partir de cela, plusieurs manipulations sont réalisées pour préparer les données aux modélisations :

  • sélection : sélection des variables qui entreront dans les modélisations ;
  • filtres : premier filtre pour ne garder que les articles pour lesquels toutes les variables à intégrer aux modèles ont une valeur connue (not NA), puis second filtre sur les articles pour lesquels des APC ont été payés (apc_has_been_paid == 1, amount_apc_EUR_trusted > 0).

Le nombre d’articles dans notre base finale pour les modélisations est ainsi de 119616. En ne gardant que les sources d’APC fiables, 28554 articles ont été écartés, soit 19.27% des articles avec des APC payés et valeurs connues pour les autres variables.


Valeurs aberrantes

graph <- bso_pop %>% ggplot(aes(x = amount_apc_EUR_trusted)) +
      geom_boxplot(fill = "#112446") +
      coord_flip() +
      labs(x = "Montant des APC fiables", y = " ", title = "Boxplot de Y") +
      theme_ipsum() +
      theme(axis.title.y = element_text(size=12), axis.title.x = element_text(size=12))
#boxplot <- ggplotly(graph)
graph1 <- bso_pop %>% ggplot(aes(x = amount_apc_EUR_trusted)) +
      geom_histogram(bins = 30L, fill = "#112446") +
      labs(x = "Montant des APC fiables", y = "Nombre d'articles", title = "Histogramme de Y") +
      theme_ipsum() +
      theme(axis.title.y = element_text(size=12), axis.title.x = element_text(size=12)) +
      stat_function(fun = dnorm, args = list(mean = mean(bso_pop$amount_apc_EUR_trusted), sd = sd(bso_pop$amount_apc_EUR_trusted)), col = "white")
#histo <- ggplotly(graph1)
grid.arrange(graph, graph1, ncol = 2, nrow = 1)

L’étendue (écart entre les valeurs minimales et maximales) du montant des APC est importante, avec des valeurs allant de 2.27 à 9848.4 euros. Nous constatons 2 valeurs extrêmes sur le boxplot, où les montants sont respectivement de 9848.4 et 9028.43 euros, nous filtrons donc la base sur les montants inférieurs à 7000 euros (le 3e APC le plus élevé est de 6848.67 euros) et regardons à nouveau les graphiques d’évolution des APC par année, selon les différentes variables qualitatives.


Graphiques d’évolution

Discipline

# Suppression des valeurs aberrantes
bso_pop <- bso_pop %>% filter(amount_apc_EUR_trusted < 7000) %>% mutate(is_french_CA_bso_wos_openalex_single_lang = as.factor(is_french_CA_bso_wos_openalex_single_lang))

# Graphique / discipline
    # table
table <- bso_pop %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>%  group_by(year, bso_classification) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()
    # plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = bso_classification, colour = bso_classification)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par discipline et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Tier

# Graphique / tier
    # table
table <- bso_pop %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>%  group_by(year, tier) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()
    # plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = tier, colour = tier)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par tier et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


French CA

# Graphique / french CA
    # table
table <- bso_pop %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>%  group_by(year, is_french_CA_bso_wos_openalex_single_lang) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% 
  rename(is_french_CA = is_french_CA_bso_wos_openalex_single_lang)
    # plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_french_CA, colour = is_french_CA)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si le CA est français") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Couleur OA

# Graphique / oa_color_finale
    # table
table <- bso_pop %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>%  group_by(year, oa_color_finale) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup()
    # plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = oa_color_finale, colour = oa_color_finale)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par couleur OA et par année") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Covered by Couperin

table <- data %>% filter(apc_has_been_paid == 1 & !is.na(amount_apc_EUR_trusted)) %>% group_by(year, is_covered_by_couperin) %>% 
  summarise(montant_moyen = mean(amount_apc_EUR_trusted, na.rm = T)) %>% ungroup() %>% mutate(is_covered_by_couperin = as.character(is_covered_by_couperin))

  # Plot
graph <- ggplot(table, aes(x = year, y = montant_moyen, group = is_covered_by_couperin, colour = is_covered_by_couperin)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  geom_point(colour="white", size = .6, pch = 21, stroke = 1.5) +
  #scale_color_manual(values = c("#3a25ff", "#82888d", "#eeeaff", "#fffd33", "#48ffbc", "#f6f6f6")) +
  labs(x = "Année", y = "APC moyens", title = "Montants moyens des APC payés par année, selon si l'article est couvert par les accords Couperin") +
  theme_classic() +
  theme(legend.position = "right")
ggplotly(graph)


Normalité de distribution

# Normalité de Y
mean <- mean(bso_pop$amount_apc_EUR_trusted, na.rm = TRUE)   
sd <- sd(bso_pop$amount_apc_EUR_trusted, na.rm = TRUE)  
ggplot(bso_pop, aes(bso_pop$amount_apc_EUR_trusted)) +
      geom_histogram(aes(y=..density..),  color="black", fill = "steelblue", alpha = 0.2) +
    geom_density( color="red", size = 1, adjust = 5) +
    stat_function(fun = dnorm, colour = "Black", size = 1, args = list(mean = mean, sd = sd))  +
    xlim(c(min(bso_pop$amount_apc_EUR_trusted)-1, max(bso_pop$amount_apc_EUR_trusted)+1))+ ggtitle("Distribution du montant des APC fiables") +
  xlab("Montant des APC fiables") + ylab("Densité") +
     theme_classic()

# Test statistique de normalité
out <- JarqueBeraTest(bso_pop$amount_apc_EUR_trusted, robust = FALSE, method = "chisq")  # Y non normal (on rejette H0)
library(formattable)
df <- data.frame(value = unlist(out))
tdf <- as.data.frame(t(df))
formattable(tdf)
statistic.X-squared parameter.df p.value method data.name
value 47076.3232661981 2 0 Jarque Bera Test bso_pop$amount_apc_EUR_trusted

Le montant des APC fiables n’est pas normalement distribué.


Corrélations et dépendances entre les variables

# Fonction pour sauvegarder les résultats des tests
save_results <- function(output, num){
    df <- data.frame(value = unlist(output))
    tdf <- as.data.frame(t(df))
    assign(glue("tdf_{num}"), tdf[,1:5], envir = .GlobalEnv)
}


# Test de corrélation (Spearman car distribution non normale)
output <- cor.test(bso_pop$amount_apc_EUR_trusted, bso_pop$year, method = "spearman") # corrélation significative
df <- data.frame(value = unlist(output))
tdf <- as.data.frame(t(df))
tdf <- tdf %>% mutate(data.name = str_replace_all(data.name, "_", "\\_"),  #escape special characters
                          data.name = gsub("bso_pop\\$", "", data.name))
formattable(tdf)
statistic.S p.value estimate.rho null.value.rho alternative method data.name
value 242612563868318 0 0.149413478963885 0 two.sided Spearman’s rank correlation rho amount_apc_EUR_trusted and year 

# Tests de dépendance
    # Année
save_results(chisq.test(bso_pop$year, bso_pop$bso_classification), "1")
save_results(chisq.test(bso_pop$year, bso_pop$tier), "2")
save_results(chisq.test(bso_pop$year, bso_pop$oa_color_finale), "3")
save_results(chisq.test(bso_pop$year, bso_pop$is_french_CA_bso_wos_openalex_single_lang), "4")
save_results(chisq.test(bso_pop$year, bso_pop$is_covered_by_couperin), "5")
    # Discipline
save_results(chisq.test(bso_pop$bso_classification, bso_pop$tier), "6")
save_results(chisq.test(bso_pop$bso_classification, bso_pop$oa_color_finale), "7")
save_results(chisq.test(bso_pop$bso_classification, bso_pop$is_french_CA_bso_wos_openalex_single_lang), "8")
save_results(chisq.test(bso_pop$bso_classification, bso_pop$is_covered_by_couperin), "9")
    # Tier
save_results(chisq.test(bso_pop$tier, bso_pop$oa_color_finale), "10")
save_results(chisq.test(bso_pop$tier, bso_pop$is_french_CA_bso_wos_openalex_single_lang), "11")
save_results(chisq.test(bso_pop$tier, bso_pop$is_covered_by_couperin), "12")
    # Couleur OA
save_results(chisq.test(bso_pop$oa_color_finale, bso_pop$is_french_CA_bso_wos_openalex_single_lang), "13")
save_results(chisq.test(bso_pop$oa_color_finale, bso_pop$is_covered_by_couperin), "14")
    # French CA
save_results(chisq.test(bso_pop$is_french_CA_bso_wos_openalex_single_lang, bso_pop$is_covered_by_couperin), "15") #toutes dépendantes

# Présentation des résultats dans une table
output <- rbind(tdf_1,tdf_2,tdf_3,tdf_4,tdf_5,tdf_6,tdf_7,tdf_8,tdf_9,tdf_10,tdf_11,tdf_12,tdf_13,tdf_14,tdf_15) %>% 
    mutate(p.value = as.numeric(p.value),
           Dependance = case_when(p.value <= 0.05 ~ "Oui",
                                  TRUE ~ "Non"))
output <- output %>% mutate(data.name = str_replace_all(data.name, "_", "\\_"),  #escape special characters
                          data.name = gsub("bso_pop\\$", "", data.name))  #remove base name
rownames(output) <- NULL  #remove rownames
kable(output, format = "html", caption = "Résultats des tests de dépendance entre les variables qualitatives") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F, 
                html_font = "sans-serif")
Résultats des tests de dépendance entre les variables qualitatives
statistic.X-squared parameter.df p.value method data.name Dependance
1843.61062820764 63 0 Pearson’s Chi-squared test year and bso_classification Oui
10563.0744012642 21 0 Pearson’s Chi-squared test year and tier Oui
89.3806015204807 7 0 Pearson’s Chi-squared test year and oa_color_finale Oui
207.610054207291 7 0 Pearson’s Chi-squared test year and is_french_CA_bso_wos_openalex_single_lang Oui
155.89053408379 7 0 Pearson’s Chi-squared test year and is_covered_by_couperin Oui
7904.05509513946 27 0 Pearson’s Chi-squared test bso_classification and tier Oui
4499.27337899457 9 0 Pearson’s Chi-squared test bso_classification and oa_color_finale Oui
223.812265757185 9 0 Pearson’s Chi-squared test bso_classification and is_french_CA_bso_wos_openalex_single_lang Oui
13108.8454392833 9 0 Pearson’s Chi-squared test bso_classification and is_covered_by_couperin Oui
11884.0898801833 3 0 Pearson’s Chi-squared test tier and oa_color_finale Oui
228.259385498635 3 0 Pearson’s Chi-squared test tier and is_french_CA_bso_wos_openalex_single_lang Oui
15863.5696643588 3 0 Pearson’s Chi-squared test tier and is_covered_by_couperin Oui
1870.74427535178 1 0 Pearson’s Chi-squared test with Yates’ continuity correction oa_color_finale and is_french_CA_bso_wos_openalex_single_lang Oui
49845.1115001406 1 0 Pearson’s Chi-squared test with Yates’ continuity correction oa_color_finale and is_covered_by_couperin Oui
1145.57365987729 1 0 Pearson’s Chi-squared test with Yates’ continuity correction is_french_CA_bso_wos_openalex_single_lang and is_covered_by_couperin Oui


Autres tests statistiques

# Vérifier : 
    #- homoscédasticité (homogénéité de la variance au sein des sous-pops)
    #- indépendance (sous-pops doivent avoir un effet aléatoire sur Y)
    #- autocorrélation (résidus et erreurs ne doivent pas être auto-correlés)


II - Modélisations Y : 2013-2020


Péparation des données aux modélisations

Pour les modélisations, nous cherchons donc à expliquer et prédire le montant des APC fiables grâce aux variables suivantes :

  • variables explicatives (indépendantes) : oa_color_finale, tier, bso_classification, is_french_CA_bso_wos_openalex_single_lang ;
  • sources de variations (sous-populations) : year, is_covered_by_couperin.
# Train test split
train_size <- floor(0.8 * nrow(bso_pop)) # 80% of the sample size

## set the seed to make your partition reproducible
set.seed(123)
train_ind <- sample(seq_len(nrow(bso_pop)), size = train_size)

train <- bso_pop[train_ind, ]
test <- bso_pop[-train_ind, ]

Ensuite, nous séparons notre base de données en 2 échantillons : un échantillon d’apprentissage (train) et un échantillon de test (test). Le premier contient 80% des données, soit 95691 articles et le second contient les 20% derniers, soit 23923 articles. Tous les modèles seront construits sur l’échantillon d’apprentissage, puis seront validés / testés sur l’échantillon test ; s’agissant de données inconnues (les modèles n’auront pas été entraînés dessus), les appliquer à cet échantillon permettra de tester leur robustesse et leur capacité prédictive.


Régression linéaire : effets fixes


  • Modèle avec seulement l’année comme variable explicative
# Reg linéaire
lm1 <- lm(amount_apc_EUR_trusted ~ year, data = train)
summary <- tidy(lm1)
kable(summary, format = "html") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F, 
                html_font = "sans-serif")
term estimate std.error statistic p.value
(Intercept) -77428.73049 2396.139753 -32.31395 0
year 39.27705 1.187785 33.06747 0

Chaque année, le montant des APC (fiables) augmente de 39 euros en moyenne.


  • Modèle avec toutes les variables explicatives
# Reg linéaire
lm2 <- lm(amount_apc_EUR_trusted ~ year + tier + oa_color_finale + is_covered_by_couperin + is_french_CA_bso_wos_openalex_single_lang, data = train)
summary <- tidy(lm2)
kable(summary, format = "html") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F, 
                html_font = "sans-serif")
term estimate std.error statistic p.value
(Intercept) -60191.26359 2256.790367 -26.67118 0
year 30.85978 1.118668 27.58619 0
tier2 -311.08861 9.693970 -32.09094 0
tier3 -440.85370 9.574222 -46.04590 0
tier4 -446.89060 10.263922 -43.53995 0
oa_color_finalehybridOA 801.69500 8.181967 97.98316 0
is_covered_by_couperin1 -96.96469 9.262517 -10.46850 0
is_french_CA_bso_wos_openalex_single_lang1 -67.89926 4.874324 -13.92998 0


Modèles mutli-niveaux : effets fixes + aléatoires


Modèles de Joël

  • Premier modèle : random intercepts par discipline

La discipline est ici la source de variations du modèle, chaque discipline dispose de son ordonnée à l’origine.

# Modèle linéaire en utilisant le package lme4
model1_spec <- linear_reg() %>% 
    set_engine("lmer") 

# une ordonnée à l'origine qui peut varier d'une discipline à l'autre
model1_fit <- model1_spec %>% 
    fit(amount_apc_EUR_trusted ~ (1 | bso_classification), data = train) 

# vue générale du modèle
summary <- tidy(model1_fit)
kable(summary, format = "html") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F, 
                html_font = "sans-serif")
effect group term estimate std.error statistic
fixed NA (Intercept) 1495.2871 105.4812 14.17586
ran_pars bso_classification sd__(Intercept) 332.7945 NA NA
ran_pars Residual sd__Observation 810.2233 NA NA 

# les variations des APC par discipline
tidy(model1_fit, effects = "ran_vals") %>%  # on calcule les valeurs des effets aléatoires
    mutate(estimate = 1777 + estimate,  # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
           ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
           ymax = estimate + 2 * `std.error`) %>% 
    arrange(estimate) %>% 
    mutate(level = as_factor(level)) %>% 
    ggplot(aes(x = level, y = estimate)) +
    geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
    geom_hline(yintercept = 1777, linetype = 2) +
    coord_flip()


  • Deuxième modèle : random intercepts par discipline et trend dépend de l’année

Par rapport au premier modèle, on ajoute l’année comme variable explicative du modèle - qui impacte donc l’angle de la pente de régression (mais celle-ci est la même par cluster, c’est-à-dire par discipline).

## Tentons un 2e modèle, dans lequel on suppose que le montant des APC augmente au cours du temps 

# on commence par transformer la variable année
train <- train %>% 
    mutate(year = year - 2013)

model2_fit <- model1_spec %>% 
    fit(amount_apc_EUR_trusted ~ year + (1 | bso_classification), data = train) # l'effet de l'année sur l'APC est supposé le même quel que soit la discipline

# vue générale du modèle
summary <- tidy(model2_fit) # l'APC moyen augmente de 34,6€ par an
kable(summary, format = "html") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F, 
                html_font = "sans-serif")
effect group term estimate std.error statistic
fixed NA (Intercept) 1311.20846 103.378135 12.68361
fixed NA year 41.67462 1.174069 35.49587
ran_pars bso_classification sd__(Intercept) 325.72589 NA NA
ran_pars Residual sd__Observation 804.94652 NA NA 

# les variations des APC par discipline - ici les valeurs sont calculées pour 2013 (year = 0)
tidy(model2_fit, effects = "ran_vals") %>%  # on calcule les valeurs des effets aléatoires
    mutate(estimate = 1565 + estimate,  # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
           ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
           ymax = estimate + 2 * `std.error`) %>% 
    arrange(estimate) %>% 
    mutate(level = as_factor(level)) %>% 
    ggplot(aes(x = level, y = estimate)) +
    geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
    geom_hline(yintercept = 1565, linetype = 2) +
    coord_flip()


  • Troisième modèle : random intercepts par discipline et trend dépend de l’année par discipline

Par rapport au deuxième modèle, on autorise la pente liée à l’année, à varier selon la discipline (more random effects).

## Dans un 3e modèle, on suppose que le montant des APC augmente au cours du temps mais potentiellement de manière différente d'une discipline à l'autre
model3_fit <- model1_spec %>% 
    fit(amount_apc_EUR_trusted ~ year + (1 + year | bso_classification), data = train) # l'effet de l'année sur l'APC est dit aléatoire (il peut varier d'une discipline à l'autre)

# vue générale du modèle
summary <- tidy(model3_fit) # l'APC moyen augmente de 29.3€ par an, mais il y a pas mal de variabilité
kable(summary, format = "html") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F, 
                html_font = "sans-serif")
effect group term estimate std.error statistic
fixed NA (Intercept) 1439.5454903 99.17308 14.515486
fixed NA year 13.7111858 10.77209 1.272844
ran_pars bso_classification sd__(Intercept) 309.8740532 NA NA
ran_pars bso_classification cor__(Intercept).year -0.0450307 NA NA
ran_pars bso_classification sd__year 32.6505586 NA NA
ran_pars Residual sd__Observation 802.6447803 NA NA 


# les variations des APC par discipline - ici les valeurs sont calculées pour 2013 (year = 0)
tidy(model3_fit, effects = "ran_vals") %>%  # on calcule les valeurs des effets aléatoires
    filter(term %in% "(Intercept)") %>% # dans un premier temps on ne regarde que l'AC moyen par discipline
    arrange(estimate) %>% 
    mutate(level = as_factor(level)) %>% 
    mutate(estimate = 1925 + estimate,  # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
           ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
           ymax = estimate + 2 * `std.error`) %>% 
    ggplot(aes(x = level, y = estimate)) +
    geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
    geom_hline(yintercept = 1925, linetype = 2) +
    coord_flip() 

tidy(model3_fit, effects = "ran_vals") %>%  # on calcule les valeurs des effets aléatoires
    filter(term %in% "year") %>% # maintenant on ne garde que les pentes
    arrange(estimate) %>% 
    mutate(level = as_factor(level)) %>% 
    mutate(estimate = 29.3 + estimate,  # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
           ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
           ymax = estimate + 2 * `std.error`) %>% 
    ggplot(aes(x = level, y = estimate)) +
    geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
    geom_hline(yintercept = 29.3, linetype = 2) +
    coord_flip() 
# dans les humanités ou en médecine, le montant des APC diminue au fil du temps tandis qu'en sciences sociales et en physique, il augmente

tidy(model3_fit, effects = "ran_vals") %>% 
    select(-`std.error`) %>% 
    pivot_wider(names_from = term, values_from = estimate)

A tibble: 10 × 5

effect group level (Intercept) year 1 ran_vals bso_classification “Biology (fond.)” 146. 56.8 2 ran_vals bso_classification “Chemistry” 391. -38.9 3 ran_vals bso_classification “Computer and informatio… -255. 20.2 4 ran_vals bso_classification”Earth, Ecology, a… 52.3 23.3 5 ran_vals bso_classification “Engineering” 47.9 -27.1 6 ran_vals bso_classification “Humanities” -342. -31.2 7 ran_vals bso_classification “Mathematics” -542. -9.75 8 ran_vals bso_classification “Medical research” 300. 17.4 9 ran_vals bso_classification “Physical sciences, Astrono… 299. -27.2 10 ran_vals bso_classification”Social sciences” -97.9 16.5


Modèles pour l’analyse

  • Premier modèle : random intercepts par année

L’année est ici la source de variations du modèle, chaque année dispose de son ordonnée à l’origine.

# Modèle linéaire en utilisant le package lme4
model1_spec <- linear_reg() %>% 
    set_engine("lmer") 

# une ordonnée à l'origine qui peut varier d'une année à l'autre
model1_fit <- model1_spec %>% 
    fit(amount_apc_EUR_trusted ~ 1 + (1 | year), data = train) 

# vue générale du modèle
summary <- tidy(model1_fit)
kable(summary, format = "html") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F, 
                html_font = "sans-serif")
effect group term estimate std.error statistic
fixed NA (Intercept) 1771.9066 41.84011 42.34947
ran_pars year sd__(Intercept) 118.0731 NA NA
ran_pars Residual sd__Observation 817.8890 NA NA 

# les variations des APC par année
tidy(model1_fit, effects = "ran_vals") %>%  # on calcule les valeurs des effets aléatoires
    mutate(estimate = 1772 + estimate,  # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
           ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
           ymax = estimate + 2 * `std.error`) %>% 
    arrange(estimate) %>% 
    mutate(level = as_factor(level)) %>% 
    ggplot(aes(x = level, y = estimate)) +
    geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
    geom_hline(yintercept = 1772, linetype = 2) +
    coord_flip()


  • Deuxième modèle : random intercepts par année et trend dépend de la discipline

Par rapport au premier modèle, on ajoute la discipline comme variable explicative du modèle - qui impacte donc l’angle de la pente de régression (mais celle-ci est la même par cluster, c’est-à-dire par année).

## Tentons un 2e modèle, dans lequel on suppose que le montant des APC dépend de la discipline
model2_fit <- model1_spec %>% 
    fit(amount_apc_EUR_trusted ~ bso_classification + (1 | year), data = train) # l'effet de la discipline sur l'APC est supposé le même quelle que soit l'année

# vue générale du modèle
summary <- tidy(model2_fit)
kable(summary, format = "html") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F, 
                html_font = "sans-serif")
effect group term estimate std.error statistic
fixed NA (Intercept) 1849.35886 43.361568 42.649723
fixed NA bso_classificationChemistry -197.16272 15.207234 -12.965061
fixed NA bso_classificationComputer and information sciences -553.73092 19.908469 -27.813838
fixed NA bso_classificationEarth, Ecology, Energy and applied biology -232.68578 9.258577 -25.131917
fixed NA bso_classificationEngineering -486.09551 27.113104 -17.928434
fixed NA bso_classificationHumanities -882.32097 34.474890 -25.593148
fixed NA bso_classificationMathematics -925.75223 34.972179 -26.471106
fixed NA bso_classificationMedical research -13.18786 6.038334 -2.184023
fixed NA bso_classificationPhysical sciences, Astronomy -227.31944 10.420008 -21.815669
fixed NA bso_classificationSocial sciences -417.68790 34.692059 -12.039871
ran_pars year sd__(Intercept) 122.01518 NA NA
ran_pars Residual sd__Observation 804.11607 NA NA 


# les variations des APC par année - ici les valeurs sont calculées pour la discipline "Social"
tidy(model2_fit, effects = "ran_vals") %>%  # on calcule les valeurs des effets aléatoires
    mutate(estimate = 1852 + estimate,  # pour rajouter l'ordonnée à l'origine pour l'ensemble de la population
           ymin = estimate - 2 * `std.error`, # pour calculer des intervalles à 95%
           ymax = estimate + 2 * `std.error`) %>% 
    arrange(estimate) %>% 
    mutate(level = as_factor(level)) %>% 
    ggplot(aes(x = level, y = estimate)) +
    geom_pointrange(aes(ymin = ymin, ymax = ymax)) +
    geom_hline(yintercept = 1852, linetype = 2) +
    coord_flip()


  • Troisième modèle : random intercepts par année et trend dépend de la discipline par année

Par rapport au deuxième modèle, on autorise la pente liée à la discipline, à varier selon l’année (more random effects).

### Hypothèses

# Test d'homogénéité de la variance (H0)
bartlett.test(amount_apc_EUR_trusted ~ interaction(year, is_covered_by_couperin), data = bso_pop)
Bartlett test of homogeneity of variances

data: amount_apc_EUR_trusted by interaction(year, is_covered_by_couperin) Bartlett’s K-squared = 1959.9, df = 15, p-value < 2.2e-16

    #-->use weights arguments from lme()

# Test d'autocorrélation des résidus

# Taille des sous-populations
table(bso_pop$is_covered_by_couperin) # échantillons déséquilibrés
 0      1

102667 16947

table(bso_pop$year) # de 8184 à 26054 (peu d'obs pour 2013-2014 car suppression des NA)

2013 2014 2015 2016 2017 2018 2019 2020 8184 9537 11251 12978 14839 16943 19828 26054

    # 16 clusters
    #--> passer en bayésien


III - Imputation des valeurs manquantes de Y : 2013-2020

Homogénéité des articles avec Y inconnus

# Filtre sur les Y is_french_CA manquants
sample <- data %>% filter(is.na(is_french_CA_bso_wos_openalex_single_lang))


IV - Prévisions Y : 2021-2030


Scénario 1


Scénario 2


Scénario 3

 

Document sous licence ouverte